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Combining  nimierical  techniques  with  ideas  from  symbolic  computation 
and  with  methods  incorporating  knowledge  of  science  and  mathematics  leads 
to  a  new  category  of  intelligent  computationzd  tools  for  scientists  and  en¬ 
gineers.  These  tools  autonomously  prepare  simulation  experiments  from 
high-level  specifications  of  physical  models.  For  computationally  intensive 
experiments,  they  automatically  design  specicd-purpose  numerical  engines 
optimized  to  perform  the  necessary  computations.  They  actively  monitor 
numerical  and  physical  experiments.  They  interpret  experimental  data  and 
formulate  numerical  results  in  qualitative  terms.  They  enable  their  human 
users  to  control  computational  experiments  in  terms  of  high-level  behavioral 
descriptions. 

As  an  example  of  such  a  tool,  imagine  an  ocean  engineer  designing  a 
offshore  mooring  tower  for  large  ships.  When  standing  free,  such  a  tower  can 
be  modeled  in  a  straightforward  way  as  an  inverted  pendulum  anchored  to 
the  sea  bed,  driven  by  wave  motion,  and  restored  to  vertical  position  by  its 
buoyancy  in  sea  water.  However,  a  meissive  ship  attached  to  the  mooring 
structure  vastly  complicates  the  dyneimics:  when  the  mooring  line  is  slack, 
the  restoring  force  is  the  same  as  before;  but  when  the  line  is  taut,  the  ship 
also  exerts  a  force  through  tension  in  the  line.  Numerical  simulations  of  such 
systems  have  shown  that  the  effect  of  even  small  waves  cein  leeid  to  potentially 
dangerous  resonances  that  can  make  the  structure  fail  [25].  The  goal  of  the 
engineer  is  to  design  a  structure  that  can  withstamd  the  worst  conditions  that 
cau  be  expected  to  occur  in  the  environment  without  making  the  structure 
excessively  expensive. 

The  following  excerpt  of  a  dialogue  shows  how  we  imagine  an  intelligent 
computer  assistant  aiding  an  ocean  engineer  in  the  design  of  such  a  structure. 
Assume  that  the  engineer  has  already  proposed  a  particular  structure,  either 
by  drawing  it  or  by  combining  previously  defined  structures. 
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Computer:  The  dynamical  equations  for  this  model  are:  (displays  a 
very  large  set  of  equations). 

Engineer:  Call  that  the  full  model.  Try  to  find  a  simpler  approxi¬ 
mate  model. 

Computer:  The  time  constants  for  the  bending  moments  and  tor¬ 
sional  stresses  in  equations  1  through  25  are  very  small  compared 
with  the  period  of  lowest  mode  of  the  structure  seen  as  a  pendulum 
and  with  the  periqds  of  the  largest  components  in  our  standard  wave 
spectrum  for  this  region  of  the  gulf  coast.  Because  the  deflections  of 
the  structure  are  small  we  can  use  the  following  approximate  model: 
(displays  a  much  smaller  number  of  equations). 

Engineer:  Call  that  the  simplified  model. 

Engineer:  Characterize  the  response  of  the  simplified  model  to 
the  standard  wave  spectrum. 

Computer:  State-space  analysis  shows  that  the  response  of  this  sys¬ 
tem  can  be  broken  into  three  re^ons.  For  very  small  deflection  angles 
(smaller  than  2  degrees)  there  are  stable  oscillations,  damped  with 
time  constants  similar  to  the  oscillation  frequencies.  Significant  de¬ 
flections  toward  the  boat  show  trajectories  with  competing  domains 
of  attraction.  There  are  large  subharmonic  responses.  For  deflections 
larger  than  about  10  degrees  we  find  complex  chaotic  motions,  but 
such  large  deflections  violate  the  assumptions  of  the  model. 

Engineer:  Characterize  the  response  to  the  dominant  frequency  in 
the  spectrum  as  the  amplitude  varies  from  1  foot  to  20  foot  seas. 

Computer:  There  are  two  families  of  periodic  motions  here,  A  and 
B.  Family  A  is  present  at  the  start  of  the  parameter  range  s  =  1, 
with  a  periodic  motion  of  order  1.  At  s  =  5.2  this  motion  undergoes 
a  symmetry-brealcing  bifurcation  and  splits  into  a  pair  of  motions  Ai 
and  A2.  As  s  increeises,  both  Ai  and  A2  undergo  period-doubling 
cascades:  order  2  (s  =  6.5),  order  4  (s  =  9.8),  order  8  (s  =  10.7), 
order  16  (s  =  10.9)  The  cascade  was  not  traced  passed  the  order  16 
motions.  Chaotic  motions  were  observed  at  s  =  12.  Family  B  begins 
with  a  motion  of  order  1  arising  from  a  fold  bifurcation  at  s  =  8.  This 
motion  persists  through  the  end  of  the  parameter  range  s  =  20. 

Engineer:  How  do  the  maximum  stresses  affect  the  structure? 


Computer:  Family  A  motions  are  within  the  safe  operating  region 
for  structures  of  type  UG/751  construction,  but  the  family  B  motions 
can  easily  exceed  the  safe  operating  region  in  a  50-yeaj  storm. 

Engineer:  Perform  the  same  analysis,  assuming  that  we  triple  the 
bouyancy  of  the  float. 

Computer:  The  qualitative  behavior  in  the  two  situations  differs 
slightly  in  the  initial  phase  of  growth  of  oscillations.  Once  oscillations 
are  established  the  behavior  is  approximately  the  same.  The  maxi¬ 
mum  deflections  differ  by  less  than  10%  between  the  two  situations, 
but  the  frequencies  are  very  different. 

Engineer:  Begin  design  of  an  active  stabilizer  to  damp  the  fam¬ 
ily  B  motions.  Assume  we  install  thrusters  at  positions  halfway  up 
the  structure  and  at  90%  of  the  height  of  the  structure.  Expand  the 
full  model  to  include  drives  from  these  thrusters.  . . . 


No  systems  as  sophisticated  as  this  engineer’s  assistant  yet  exist.  How¬ 
ever,  computer  programs  incorporating  important  pieces  of  it  are  already  in 
place.  Generating  stylized  natural  language  is  not  too  hard,  although  under¬ 
standing  unstructured  English  text  is  still  difficult,  and  we  will  not  address 
that  problem  in  this  paper.  Additionally,  our  discussion  is  not  really  about 
ocean  engineering;  the  scenario  above  is  not  intended  to  illustrate  good  design 
practice  in  that  domain.  Rather,  our  concern  here  is  with  the  development  of 
intelligent  techniques  appropriate  for  the  automatic  preparation,  execution, 
and  control  of  numerical  experiments,  and  with  the  automatic  interpretation 
of  their  results. 

•  Our  envisioned  engineer’s  assistant  begins  with  a  description  of  a  mech¬ 
anism  and  automatically  generates  efficient  numeric2d  programs  that 
predict  its  dynamical  behavior.  This  may  require  more  tham  just  straight¬ 
forward  simulation.  A  stability- amalysis  task  such  m  “characterize  the 
response  to  the  dominant  frequency  in  the  spectrum”  requires  com¬ 
piling  procedures  that  evolve,  in  addition  to  the  state,  the  variations 
with  respect  to  changes  in  initial  conditions  and  the  sensitivities  with 
respect  to  changes  in  parameters. 
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•  The  engineer’s  assistant  automatically  prepares  high-performzmce  nu¬ 
merical  experiments.  It  has  extensive  knowledge  of  numerical  meth¬ 
ods  and  it  can  compose  appropriate  and  correct  numerical  procedures 
tailored  to  the  specific  application.  For  critical  applications,  this  com¬ 
pilation  can  be  targeted  to  the  automatic  synthesis  of  special-purpose 
hardware. 

•  The  engineer’s  assistant  interprets  the  results  of  numerical  experiments 
in  high-level  qualitative  terms.  This  interpretation  is  based  on  general 
mathematical  and  physical  knowledge  that  constrains  the  kind  of  be¬ 
havior  to  expect.  The  interpretation  is  used  to  prepare  a  report  to  the 
user,  but  it  is  also  used  in  the  experimental  protocol.  The  surmnary  of 
behavior  produced  from  observations  of  the  results  of  previous  exper¬ 
iments  is  used  to  automatically  select  critical  values  of  experimental 
parameters  for  subsequent  experiments,  thus  efficiently  uncovering  the 
salient  phenomena. 

Section  one  of  this  paper  demonstrates  significant  portions  of  these  ca¬ 
pabilities.  These  include  the  automatic  preparation  and  monitoring  of  nu¬ 
merical  simulations,  the  automatic  generation  of  qualitative  interpretations 
of  numerical  results,  and  the  achievement  of  breakthrough  performance  on 
computationally-demanding  problems  with  the  aid  of  specially-designed  com¬ 
puters.  (Our  special-purpose  engine  for  computing  planetary  motions  has 
produced  the  first  solid  numerical  evidence  that  the  solar-system’s  long-term 
dynamics  is  cheuitic,  thereby  answering  the  famous  question  of  the  stability 
of  the  solar  system.) 

Section  two  takes  a  closer  look  at  the  technology  behind  these  demonstra¬ 
tion  results.  We  explain  how  algorithms  from  computer  vision  are  applied  to 
interpret  phase-space  diagrams  in  dynamics.  We  illustrate  how  knowledge 
about  dynamical  systems  can  be  encoded  using  constraints  and  symbolic 
rules.  We  show  how  to  formulate  numerical  algorithms  at  appropriate  levels 
of  abstraction  with  higher-order  procedures  and  how  to  combine  these  with 
symbolic  algebra  to  automatically  generate  numerical  programs. 

Section  three  sketches  some  next  steps  required  to  realize  the  vision  of 
systems  like  the  engineer’s  assistant. 


1  Numerical  modeling  can  be  automated 


In  a  t3rpicaJ  numerical  modeling  study,  an  investigator  repeatedly  prepares 
and  runs  a  series  of  computations  and  examines  the  results  at  each  step  to 
select  interesting  new  values  for  parameters  and  initiad  conditions.  When 
enough  values  have  been  tried,  the  investigator  classifies  and  interprets  the 
results.  Even  with  powerful  numerical  computers,  this  process  requires  sub¬ 
stantial  human  effort  to  prepare  simulations,  eind  it  rehes  upon  significant 
human  judgment  to  choose  interesting  values  for  parameters,  to  determine 
when  a  simulation  run  is  complete,  amd  to  interpret  numerical  results  in 
qualitative  terms. 

This  section  exhibits  three  programs  that  automate  much  of  the  above 
process.  The  Bifurcation  Interpreter  investigates  the  steady-state  orbits  in 
parameterized  families  of  dynamical  systems,  classifying  the  types  of  orbits 
and  the  bifurcations  through  which  they  change  as  parameters  vary.  The 
KAM  program  autonomously  explores  nonlinear  conservative  systems  and 
produces  qualitative  descriptions  of  phase-space  portraits  and  bifurcations. 
Both  programs  automatically  generate  summary  reports  similar  to  those  ap¬ 
pearing  in  published  papers  in  the  experimental  dynamics  literature  and  in 
engineering  studies  of  artifacts  that  have  complex  dynamics,  such  as  mrfoils, 
ship  hulls,  and  mooring  structures.  In  addition,  the  capabilities  demonstrated 
by  these  programs  have  application  in  the  design  of  intelligent  automatic  con¬ 
trol  systems.  The  breadth  of  applicability  is  illustrated  by  the  Kineticist’s 
Workbench,  a  program  that  models  how  chemists  understand  complex  chem¬ 
ical  reactions.  It  combines  numerical  and  symbolic  methods  to  characterize 
reaction  mechanisms  in  qualitative  terms  that  axe  useful  for  the  working 
chemist. 

We  also  discuss  the  plaice  of  special-purpose  numericad  engines  as  scientific 
instruments  amd  survey  significant  results  in  planetary  dynamics  obtained 
using  the  Digital  Orrery. 
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1.1  Programs  can  discover  and  interpret  qualitative 
behavior 

In  a  nonlinear  dynamical  system  with  a  periodic  drive,  motion  starting  from 
any  set  of  initial  conditions  will  typically  evolve  to  a  steady-state  orbit.  ^ 
For  a  parameterized  family  of  dynamical  systems,  tracing  the  changes  in 
steady-state  orbits  as  the  parameters  vary  provides  a  valuable  sununaxy  of 
the  family’s  quahtative  behavior.  Much  research  in  nonlinear  dyniimics  is 
devoted  to  studying  these  bifurcations,  or  changes  in  type,  of  steady-state 
orbits.  For  one-parameter  families  at  least,  the  bifurcations  genericaJly  en¬ 
countered  have  been  classified  and  are  well- understood.  Some  examples  are 
the  fold  bifurcation,  at  which  a  stable  orbit  can  appear  or  vanish,  the  flip 
bifurcation,  at  which  the  period  of  an  orbit  doubles,  and  the  pitchfork  bifur¬ 
cation,  at  which  an  orbit  splits  into  two  orbits  of  the  same  period.  There  are 
also  conunonly- observed  bifurcation  sequences  that  occur  as  the  parzuneter 
vMies.  An  example  is  the  period- doubling  cascade,  where  the  order  of  an 
orbit  successively  doubles  via  a  sequence  of  increasingly  closely  spaced  flip 
bifurcations,  producing  chaos.^ 

Dynamicists  commonly  gain  insight  into  the  qualitative  behavior  of  non¬ 
linear  systems  by  developing  summary  descriptions  of  steady-state  orbits 
and  bifurcations.  Figure  1  reproduced  here  from  [7]  shows  a  schematic  sum¬ 
mary  drawn  by  a  physicist  based  on  numerical  studies  of  the  two-dimensional 
Navier-Stokes  equation  for  an  incompressible  fluid.  As  the  Reynolds  number 
of  the  fluid  increases,  the  steady-state  orbits  evolve  through  a  sequence  of  bi¬ 
furcations.  The  diagreun  summarizes  how  the  evolving  orbits  can  be  grouped 
into  four  distinct  families. 

The  Bifurcation  Interpreter,  a  computer  program  being  developed  at  MIT 
by  II.  Abelson,  autoiaaticaliy  generates  such  summary  descriptions  for  one- 
parameter  families  of  periodically-driven  dynamical  systems.  The  dynamical 

^Possible  types  of  steady-state  orbits  are  periodic  orbits,  quasi-periodic  orbits  (which 
have  discrete-frequency  spectra,  but  not  at  rational  multiples  of  the  drive  period),  and 
chaotic  orbits  (which,  loosely  speaking,  are  steady-state  orbits  that  are  neither  periodic 
not  quasi-periodic). 

^Various  authors  use  different,  and  sometimes  incompatible  terminology  to  refer  to 
these  bifurcation  types.  For  example,  the  flip  is  sometimes  called  a  cusp  or  a  pitchfork. 
We  have  adopted  the  terminology  used  in  the  book  by  Thompson  and  Stewart  [26],  which 
provides  an  introduction  to  the  methods  of  nonlinear  dynamics  together  with  an  extensive 
bibliography. 
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Figure  1;  This  diagram,  reproduced  from  a  published  paper  in  fluid  mechanics,  is  a 
physicist’s  schematic  summary  description  of  an  approximation  to  the  two-dimensional 
Navier-Stokes  equation  for  an  incompressible  fluid.  The  varying  parameter  here  is  the 
Reynolds  number. 


Figure  2;  The  Bifurcali^^n  Interpreter  automatically  generates  summary  descriptions 
of  dynamical  systems  similar  to  those  appearing  in  published  papers.  Ihis  (preliminary 
version)  diagram  describes  the  behavior  of  a  periodically-impacted  hinged  beam  as  the 
load  varies,  exhibiting  the  evolution  of  two  different  families  of  steady-state  motions. 
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"  *  slem  can  be  specified  by  differential  equations  to  be  integrated,  by  a  pe¬ 
riod  map  that  directly  computes  successive  states  at  multiples  of  the  drive 
period,  or  by  a  description  of  a  physical  model  such  cis  an  electrical  network. 
Given  a  dynamical  system,  a  parameter  range  to  explore,  and  a  domain  in 
state-space,  the  interpreter  discovers  periodic  orbits,  tracks  their  evolution 
as  the  parameter  varies,  and  locates  and  classifies  bifurcations.  Using  this 
information,  the  program  categorizes  the  orbits  into  families  and  produces  a 
summary  report  that  describes  each  family  and  its  evolution  through  bifur¬ 
cations. 

Here  is  a  sample  input  to  the  interpreter: 


Period  map 


1  -  ■  1  -  e' 

^  2m  '  y 

y  ~e  ^”*asinx -f 


Values  for  fixed  parameters 


m 


O.Itt 


Range  for  varying  parameter  a  1  to  7 


Bounds  on  state  variables 


X 

— TT  to  TT  (periodic) 

y 

—4  to  4 

The  system  to  be  investigated  models  the  vibration  of  a  hinged  bar  with 
viscous  damping  subjected  to  a  fixed-direction  periodic  impact  load  at  the 
free  end.  This  problem  is  discussed  in  [13],  which  derives  the  period  map  for 
this  motion.  Here  i  is  the  angular  displacement  and  y  is  the  angular  velocity, 
a  specifies  magnitude  of  the  load  and  m  specifies  the  damping  factor.  The 
interpreter  is  asked  to  explore  the  system  as  the  load  ranges  from  1  to  7. 
With  this  specification,  the  interpreter  analyzes  the  system  and  generates 
the  following  report: 

There  are  2  distinct  families  of  periodic  orbits,  A  and  B. 

Family  A  is  already  present  at  the  start  of  the  parameter  range  a  =  1 
as  a  periodic  orbit  Ao  of  order  1.  \X  a  —  4.130  there  is  a  supercritical 
flip  bifurcation  at  which  A©  undergoes  period  doubling  to  produce 


a  periodic  orbit  Ai  of  order  2.  At  a  =  6.489  there  is  a  supercritical 
pitchfork  bifurcation  at  which  the  family  A  splits  into  subfamilies  A(l) 
and  i4(2),  beginning  with  Ai  splitting  into  two  periodic  orbits  of  order 
2.  As  the  parameter  a  increases,  each  subfamily  undergoes  a  period¬ 
doubling  cascade  via  a  sequence  of  supercritical  flip  bifurcations  to 
order  4  at  a  =  6.838,  order  8  at  a  =  6.891,  order  16  at  a  =  6.901, 
order  32  at  a  =  6.903.  The  period- doubling  cascade  was  not  traced 
past  the  order  32  orbit,  which  apparently  period  doubles  again  at 
a  =  6.904. 

Family  B  first  appears  at  a  =  3.969  with  an  orbit  Bo  of  order  4 
appearing  at  a  fold  bifurcation.  As  the  parameter  a  increases,  B 
undergoes  a  period-doubling  cascade  via  a  sequence  of  supercritical 
flip  bifurcations  to  order  8  at  a  =  4.239,  order  16  at  o  =  4.239,  order 
32  at  a  =  4.251.  The  period-doubling  cascade  was  not  traced  past  the 
order  32  orbit,  which  apparently  period  doubles  again  at  a  =  4.252. 

The  program  can  display  this  same  information  as  a  diagram  (figure  2) 
that  is  similar  in  style  to  manually- developed  Navier-Stokes  analysis  in  fig¬ 
ure  1.^ 


1.2  Smart  programs  can  see  what  not  to  compute 

Dynamical  behavior  is  complex,  but  it  is  not  arbitrary.  There  is  structure 
on  phase  space  that  restricts  the  classes  of  legal  trajectories  and  provides  a 
grammar  of  legal  phase  portraits.  For  example,  trajectories  of  autonomous 
systems  cannot  intersect,  and  an  we  vary  the  initial  conditions,  the  trajec¬ 
tories  vary  smoothly  except  at  isolated  places  where  the  behavior  changes. 
As  we  vary  parameters,  the  phase  portrait  changes  qualitatively  only  at  bi¬ 
furcations.  In  Hamiltonian  systems  the  evolution  of  the  pheise  space  is  area- 
preserving,  which  greatly  restricts  the  cl£isses  of  possible  structures  that  can 
occur  in  the  phaise  space.  This  kind  of  knowledge  enables  dynamicists  to  infer 
a  good  understanding  of  a  physical  system  from  only  a  small,  but  well-chosen, 
set  of  experiments. 

The  phase  portrait  in  figure  3,  taken  from  a  historically  important  pa¬ 
per  in  dynamics  by  M.  Henon  [9],  describes  how  adding  a  simple  quadratic 

^The  diagram-generation  program  illustrated  in  figure  2  was  developed  by  Ognen 
Nastov. 
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X  I— »•  xcosa  —  (y  —  x^)sina 
y  xsina  +  (y  —  x^)  cosa 


Figure  3:  This  phase  portrait  is  Henon’s  summary  of  the  dynamics  of  the  map  for 
cos  a  =  .24. 

nonlinearity  to  a  linear  rotation  can  lead  to  dramatic  changes  in  dynamical 
behavior.  Observe  how  the  figure  characterizes  the  dynamics  by  showing 
only  a  few  orbits.  Presumably,  Henon  was  able  to  generate  this  figure  after 
performing  only  a  few  judiciously  chosen  numeric£il  experiments. 

The  KAM  program  developed  by  K.  Yip  at  MIT  can  analyze  systems  in 
the  same  way  [29,  30].  It  knows  enough  about  the  constraints  on  the  structure 
of  phase  space  to  choose  initial  conditions  and  parameters  as  cleverly  as  an 
expert  dynamicist.  RAM’s  summary  description  of  Henon’s  map  is  shown 
in  figure  4.  Observe  that  this  is  almost  identical  to  the  summary  presented 
by  Henon.  Moreover,  RAM  was  able  to  deduce  this  description  after  trying 
only  ten  initial  conditions. 

RAM’s  ability  to  control  numerical  experiments  arises  from  the  fact  that 
it  not  only  produces  pictures  for  us  to  see — it  also  looks  at  the  pictures 
it  draws,  visually  recognizing  and  classifying  different  orbit  types  as  they 
numerically  evolve.  By  combining  techniques  from  computer  vision  with 
sophisticated  dynamical  invariants,  RAM  is  able  to  exploit  mathematical 
knowledge,  represented  in  terms  of  a  “grammar”  that  dictates  consistency 
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Figure  4;  The  KAM  program  generates  this  summary  picture  of  Henon’s  map. 

constraints  on  the  structure  of  phase  space.  When  it  chooses  new  initial 
conditions  to  explore,  it  does  so  in  an  attempt  to  make  the  picture  consistent 
with  these  constraints.  In  addition  to  drawing  the  picture,  KAM  generates 
a  textual  analysis  that  explains  what  the  program  “sees.”  Here  is  KAM’s 
description  of  the  picture  it  generates  for  Henon’s  map. 

The  portrait  has  an  elliptic  fixed  point  at  (0,0).  Surrounding  the 
fixed  point  is  a  regular  region  bounded  by  a  KAM  curve  with  rotation 
number  between  1/5  and  1/4.  Outside  the  regular  region  lies  a  chain  of 
5  islands.  The  island  chain  is  bounded  by  a  KAM  curve  with  rotation 
number  between  4/21  and  5/26.  The  outermost  region  is  occupied  by 
chaotic  orbits  that  eventually  escape. 


1.3  Programs  can  construct  and  analyze  approxima¬ 
tions 

A  powerful  strategy  for  analyzing  a  complicated  dynamical  system  is  to  ap¬ 
proximate  it  with  a  simpler  system,  analyze  the  approximation,  and  map 
the  results  back  to  the  original  system.  The  approximations  must  be  accu¬ 
rate  enough  to  reproduce  the  essential  properties  of  the  original  system,  yet 
simple  enough  to  be  analyzed  efficiently.  Human  experts  have  found  that 
piecewise  linear  approximations  satisfy  both  criteria  for  a  wide  class  of  mod¬ 
els.  The  PLR  program,  developed  by  E.  Sacks  at  MIT,  exploits  this  fact  to 
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automate  the  analysis  of  second-order  autonomous  ordinary  differential  equa¬ 
tions  [21,  22].  It  derives  the  qualitative  behavior  of  intractable  equations  by 
approximating  them  with  piecewise  linear  equations  and  constructing  phase 
diagrams  of  the  approximations. 

Plr  constructs  a  composite  phase  diagram  for  a  piecewise-linear  system 
by  combining  the  local  phase  diagrams  of  its  linear  regions.  It  employs  the 
standard  theory  of  linear  equations  to  ascertain  the  local  phaise  diagrams. 
Linear  systems  have  simple  well-understood  dynamics.  Either  all  trajectories 
are  periodic,  all  approach  a  fixed  point,  or  all  approach  infinity.  Plr  pastes 
together  the  local  phase  diagrams  by  determining  which  sequences  of  regions 
trajectories  can  traverse.  It  summarizes  the  results  by  a  transition  graph 
whose  nodes  and  links  represent  regions  and  transitions.  Each  path  through 
the  transition  graph  of  a  piecewise- linear  system  indicates  that  trajectories 
traverse  the  corresponding  regions  in  the  prescribed  order.  Loops  denote 
trajectories  that  remain  in  one  region  forever,  whereas  longer  cycles  denote 
trajectories  that  continually  shift  between  a  sequence  of  regions. 

As  a  simple  example,  Plr  can  qualitatively  analyze  the  behavior  of  an 
undriven  van  der  Pol  oscillator,  a  simple  nonlinear  circuit  consisting  of  a 
capacitor,  an  inductor,  and  a  nonlinear  resistor  connected  in  series.  The 
current  through  the  circuit  obeys  the  equation 

+  +  =  (1) 

with  C  the  capacitance,  L  the  inductance,  and  k  a  scaling  factor.  Plr 
approximates  this  equation  with  a  piecewise  linear  equation  and  constructs 
the  phase  diagram  and  transition  graph  shown  in  Figure  5.  It  deduces  that 
the  system  oscillates  from  the  fact  that  tracing  edges  starting  from  any  node 
in  the  graph  leads  to  a  cycle.  Intuitively,  the  system  oscillates  because  the 
nonlinear  resistor  adds  energy  to  the  circuit  at  low  currents  and  drains  energy 
at  high  currents. 

1.4  Domain  knowledge  can  guide  numerical  modeling 

M.  Eisenberg’s  Kineticist’s  Workbench,  also  being  developed  at  MIT,  is  a 
program  that  combines  general  knowledge  of  dynamics  with  specific  knowl¬ 
edge  about  chemical  reactions  in  the  analysis,  understanding,  and  simulation 
of  complex  chemical  reaction  mechanisms. 


Figure  5:  Plr’s  phase  diagram  and  transition  graph  for  its  piecewise  linear 
van  der  Pol  approximation.  Arrows  indicate  boundaries  that  trajectories 
cross. 

Chemists,  in  trying  to  model  reactions,  typically  hypothesize  a  set  of  ele¬ 
mentary  reaction  steps  {corresponding  to  molecular  collisions)  that  constitute 
a  proposed  pathway  for  the  overall  reaction.  This  collection  of  elementary 
steps  may  be  large.  It  usually  gives  rise  to  a  mathematical  model  consist¬ 
ing  of  many  tightly- coupled  nonlinear  diflFerential  equations.  The  problem  of 
simulating  such  a  system  can  be  formidable,  but  a  simulation  merely  pro¬ 
vides  numerical  results.  Even  more  important  to  the  chemist  is  to  achieve 
some  sort  of  qualitative  understanding  of  the  reaction  mechanism.  The  Ki- 
neticist’s  Workbench  combines  numerical  simulation  with  a  variety  of  sym¬ 
bolic  methods  in  order  to  characterize  reaction  mechanisms  in  terms  that  are 
meaningful  to  the  working  chemist  such  as  “equilibrium,”  “steady  states,” 
and  “oscillations.” 

The  Workbench  also  suggests  ways  to  simplify  complex  mechanisms.  The 
key  idea  here  is  that  the  program  can  use  its  qualitative  characterization  of  a 
given  mechanism’s  behavior  to  determine  which  steps  within  the  mechanism 
have  an  important  behavioral  effect,  and  which  are  superfluous  and  may 
therefore  be  dropped  from  the  chemist’s  model. 

The  following  scenario  shows  how  the  Workbench  is  intended  to  interact 
with  a  user  in  the  analysis  of  a  chemical  process.  Tiie  initial  information  is 
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specified  as  follows: 


Step  no. 

Step 

Rate  Constant 

(1) 

C  B 

300 

(2) 

B-^C 

100 

(3) 

A-^  X 

1 

(4) 

B->rX  -^Y  YD 

1 

(5) 

2X-^Y  -*ZX 

0.333 

(6) 

X  E 

1 

(7) 

C  F 

100 

External  sources  for  species  C  (rate  of  addition;  100  mmol/liter-sec) 

Concentration  of  A  assumed  constant  at  1  mmol/liter.  Initial  concen¬ 
trations  of  all  species  besides  A  are  0  mmol/liter. 


With  this  specification,  the  Workbench  responds: 


The  mechanism  may  be  decomposed  into  two  submechanisms, 
first  may  be  written  as  follows:"* 


The 


(1.1)  [External]  ^  C  B -*  [External] 


This  submechanism  results  in  a  rapid  equilibrium  between  C  and  B\ 
the  final  concentrations  of  both  C  and  5(1  and  3  mmol/liter,  respec¬ 
tively)  are  independent  of  their  initial  concentrations.® 

The  second  submechanism  may  be  written  as  follows: 

(2.1)  [External]  X 

(2.2)  BAX^Y 

(2.3)  2X-\-Y  -*ZX 

■'The  notation  [Exiemal]  used  in  this  example  is  a  catchall  used  to  indicate  the  presence 
of  external  sources  or  sinks  for  some  species.  The  notation  corresponds  roughly  to  the  “zero 
complex”  notation  of  Feinberg  [6]. 

®In  performing  this  simplification,  the  Workbench  has  treated  the  original  reaction  (7) 
as  the  presence  of  an  external  sink  for  species  C,  since  the  concentration  of  species  F 
merely  grows  with  time  and  does  not  affect  any  other  species.  Similarly,  reaction  (4), 
being  slow  compared  to  all  the  reactions  included  in  (11),  may  be  treated  as  a  slow  (but 
time- varying)  sink  for  species  B.  This  reasoning  could  be  exposed  at  the  user’s  request. 


Since  the  concentration  of  B  in  the  fi.  *  submechanism  is  constant  at 
3  mmol/liter,  the  second  submechanism  .nay  be  further  simplified  as 
follows: 


(2.1')  [External]  «-*  A'  —  Y 
(2.2')  2A'  +  r  -  3A' 


Simulation  of  this  second  mechanism  reveals  that  the  concentrations  of 
X  and  Y  oscillate  over  time.  Each  oscillation  consists  of  five  phases: 

(1)  slow  growth  in  the  concentrations  of  AT  and  F;  (2)  more  rapid 
growth  in  AT;  (3)  large  jump  in  X  and  decline  in  Y;  (4)  large  decline 
in  X  and  increase  in  Y;  (5)  slower  decline  in  X  and  slow  increase  in 
Y .  The  time  between  oscillation  peaks  is  approximately  17.8  seconds. 

The  chemist  now  asks  the  program  to  compare  the  results  of  the  approxi¬ 
mate  mechanism  (2.1'-2.2')  with  the  numerically  exact  results  of  mechanism 
(1-7),  comparing  the  qualitative  behavior  of  the  two  mechanisms  as  reflected 
in  species  X  and  Y.  The  workbench  responds: 

The  qualitative  behavior  of  the  two  mechanisms  differs  slightly  during 
the  initial  phase  of  slow  growth  in  X  and  Y .  Once  oscillations  begin, 
the  exact  mechanism  (1-7)  displays  the  same  qualitative  structure  as 
the  approximate  mechanism  (2.1'-2.2'),  and  the  period  is  still  approx¬ 
imately  17.8  seconds.  The  concentrations  of  X  and  Y  at  comparable 
points  in  the  oscillations  differ  by  less  than  2  percent  between  the  two 
mechanisms. 

The  chemist  now  tries  a  different  variation  on  the  original  reaction,  and 
alters  the  original  mechanism  so  that  the  concentration  of  A  is  not  constant 
over  time,  but  is  rather  an  explicit  function  of  time,  as  given  by  the  equation 
[A]  =  1  -|-  sint.  The  workbench  responds: 

The  decomposition  into  submechanisms  is  not  affected;  nor  is  the  con¬ 
stant  concentration  of  B.  Again,  the  concentrations  of  X  and  Y  ap¬ 
pear  to  be  oscillating,  but  each  oscillation  now  consists  of  eight  phases, 
instead  of  five,  as  before. 

It  is  worth  noting  some  of  the  key  features  of  this  example.  First,  the 
Workbench  is  able  to  decompose  the  original  mechanism  into  two  submech¬ 
anisms,  each  of  which  is  capable  of  independent  simulation;  this  simplifies 
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both  the  analysis  and  simulation  of  the  larger  mechanism.  Second,  the  Work¬ 
bench  is  able  to  decompose  the  first  of  the  two  submechanisms  in  terms  of  a 
dichotomy  between  fast  and  slow  steps;  this  allows  the  program  to  approxi¬ 
mate  the  submechanism  as  a  system  in  equilibrium.  Third,  the  program  uses 
numerical  simulation  to  derive  equilibrium  concentrations  for  this  submech- 
ajiism.  Finally,  the  Workbench  is  able  to  describe  the  results  of  simulating 
the  second  submechanism  in  terms  of  a  succession  of  qualitative  episodes 
characterized  by  changing  growth  rates  of  the  species  X  and  Y . 

1.5  Fast  computers  need  not  be  large  or  expensive 

Numerical  modeling  often  requires  substantial  resources.  Scientists  and  en¬ 
gineers  have  traditionally  obtained  these  resources  either  by  acquiring  large- 
scale  computers  or  renting  time  on  them.  However,  a  specialized  computer 
can  be  simple  and  physically  small.  Indeed,  it  may  be  just  as  easy  to  design, 
build,  and  program  a  special-purpose  computer  than  to  develop  software  for 
general-purpose  supercomputers.  Moreover,  the  specialized  computer  can  be¬ 
come  an  ordinary  experimental  instrument  belonging  to  the  research  group 
that  made  it,  thus  avoiding  the  administrative  burden  and  the  scheduling 
problems  associated  with  expensive,  shared  resources. 

The  question  of  the  stability  of  the  solar  system  is  probably  the  most 
famous  longstanding  problem  in  astrodynamics.  In  fact,  it  was  investigations 
into  precisely  this  problem  that  inspired  Poincare  to  develop  the  modern 
qualitative  theory  of  dynamical  systems.  In  1988,  G.  Sussman  and  J.  Wisdom 
completed  a  series  of  numerical  experiments  at  MIT  demonstrating  that  the 
long-term  motion  of  the  planet  Pluto,  and  by  implication  the  dynamics  of 
the  Solar  System,  is  chaotic  [24]. 

The  stability  question  was  settled  using  the  Digital  Orrery  [4],  a  special- 
purpose  numerical  engine  optimized  for  high-precision  numerical  integrations 
of  the  equations  of  motion  of  small  numbers  of  gravitationally  interacting 
bodies.  Using  1980  technology,  the  device  is  about  1  cubic  foot  of  electronics, 
dissipating  150  watts.  On  the  problem  it  was  designed  to  solve,  it  is  measured 
to  be  60  times  faster  than  a  VAX  11/780  with  FPA,  or  1/3  the  speed  of  a 
Cray  1. 

Figure  6  shows  the  exponential  divergence  of  nearby  Pluto  trajectories 
over  400  million  years.  This  data  is  taken  from  an  845-million-year  integra¬ 
tion  performed  with  the  Orrery.  Before  the  Orrery,  high-precision  integra- 
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tions  over  simulated  times  of  millions  of  years  were  prohibitively  expensive. 
The  longest  previous  integration  of  the  outer  planets  was  for  five  million 
years,  performed  on  a  Japanese  supercomputer  in  1984  [14].  Even  though 
the  Orrery  is  not  as  fast  as  the  fastest  supercomputer,  its  small  scale  and 
relative  low  cost  mean  that  it  can  be  dedicated  to  long  computations  in  ways 
that  a  conventional  supercomputer  could  not.  To  perform  the  integration 
that  established  Pluto’s  chaotic  behavior,  the  Orrery  ran  continually  for  five 
months. 

The  Orrery  was  designed  and  built  by  six  people  in  only  nine  months. 
This  was  possible  only  because  of  novel  software  support  for  the  design  pro¬ 
cess.  The  simulator  for  the  Orrery  is  partially  symbolic — simulated  registers 
hold  symbolic  values  and  simulated  arithmetic  parts  combine  these  to  pro¬ 
duce  algebraic  expressions  (in  addition  to  checking  timing  and  electrical  con¬ 
straints).  This  means  that  a  successful  simulation  yields  a  simulated  memory 
containing  algebraic  expressions  that  can  be  checked  for  correctness. 
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In  d  [AU]  (variational) 
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Figure  6;  The  exponential  divergence  of  nearby  trajectories  is  indicated  by  the  average 
linear  growth  of  the  logarithms  of  the  distance  measures  as  a  function  of  time.  In  the  upper 
trace  we  see  the  growth  of  the  variational  distance  around  a  reference  trajectory.  In  the 
lower  trace  we  see  how  two  Plutos  diverge  with  time.  The  distance  saturates  near  45 AU, 
note  that  the  semimajor  axis  of  Pluto’s  orbit  is  about  40AU.  The  variational  method  of 
studying  neighboring  trajectories  does  not  have  the  problem  of  saturation.  Note  that 
the  two  methods  are  in  excellent  agreement  until  the  two-trajectory  method  has  nearly 
saturated. 
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In  d  [AU]  (two-particle) 


2  Intelligent  numerical  computing  rests  on  A1  tech¬ 
nology 

The  illustrations  above  achieve  their  impressive  results  by  bringing  symbolic 
methods  to  bear  on  the  problems  of  numerical  computation.  Some  of  these 
techniques  are  traditional  AI  methods,  which  achieve  new  power  when  they 
are  combined  with  deep  knowledge  of  dynamical  systems.  The  KAM  pro¬ 
gram,  for  example,  uses  techniques  from  machine  vision  to  recognize  and 
classify  the  relevant  geometrical  properties  of  the  trajectories.  The  Bifurca¬ 
tion  Interpreter  uses  algebraic  manipulation  and  knowledge  about  the  local 
geometry  of  bifurcations  to  automatically  generate  numerical  procedures  that 
track  periodic  orbits.  The  key  to  automatically  generating  high-performance 
numerical  algorithms  is  to  express  knowledge  of  numerical  analysis  at  an  ap¬ 
propriate  level  of  abstraction.  This  is  supported  by  a  library  of  niunerical 
methods  that  is  organized  around  the  hberal  use  of  higher-order  procedural 
abstractions.  With  this  organization,  one  constructs  sophisticated  numerical 
methods  by  mixing  and  matching  standard  components  in  well-understood 
ways.  The  resulting  programs  are  both  more  perspicuous  and  more  robust 
than  conventional  numerical  methods.  For  example,  a  procedure  by  G.  Roy- 
lance  that  automatically  generates  special  functions  has  constucted  a  Bessel- 
function  routine  that  is  40  times  more  accurate  than  the  National  Bureau  of 
Standards  approximation,  for  the  same  amount  of  computation. 

2.1  The  KAM  program  exploits  techniques  from  com¬ 
puter  vision 

Yip’s  KAM  program  is  notable  because  it  applies  judgement,  similar  to  that 
of  an  expert  dynamicist,  in  directing  the  course  of  its  numerical  experiments. 
In  making  judicious  choices  of  what  to  try  next,  KAM  must  interpret  what 
it  sees.  This  process  occurs  in  three  phases:  aggregation,  clustering,  and 
classification.  The  images  of  an  initial  point  produced  by  iterating  the  map 
forms  a  set  of  isolated  points.  This  orbit  must  be  classified.  In  Hamiltonian 
systems  there  are  three  types  of  orbits  to  distinguish.  In  a  surface  of  section, 
periodic  orbits  appear  as  isolated  points,  quasiperiodic  orbits  appear  as  closed 
curves  or  island  chains,  and  chaotic  orbits  appear  to  take  up  regions  of  2- 
dimensional  space.  KAM  must  also  aggregate  the  components  of  an  orbit  so 
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that  it  can  be  further  classified.  It  must  be  able  to  determine  the  number 
of  islands  in  an  island  chain,  as  this  number  gives  the  period  of  the  enclosed 
periodic  point.  KAM  must  be  able  to  estimate  the  centroid  and  area  enclosed 
by  a  curve  and  to  recognize  the  shape  of  a  curve.  KAM  implements  these 
abilities  with  techniques  from  computational  geometry  and  computer  vision. 

KAM  classifies  orbits  using  methods  based  on  the  Euclidean  minimal 
spanning  tree — the  tree  that  interconnects  all  the  points  with  minimal  total 
edge  length — which  it  constructs  by  means  of  the  Prim-Dijkstra  algorithm  [5]. 
For  each  sub-tree  of  the  spanning  tree,  KAM  examines  the  degree  of  each 
of  its  nodes,  where  the  degree  of  a  node  is  the  number  of  nodes  connected 
to  it  in  the  sub-tree.  For  a  smooth  curve,  the  spanning  tree  consists  of  two 
terminal  nodes  of  degree  one  and  other  nodes  of  degree  two.  For  a  point 
set  that  fills  an  area,  its  corresponding  spanning  tree  consists  of  many  nodes 
having  degree  three  or  higher  (figure  7). 

To  aggregate  points,  KAM  deletes  from  the  tree  edges  .that  are  signif¬ 
icantly  longer  than  nearby  edges,  following  an  aggregation  algorithm  sug¬ 
gested  by  Zahn  [31].  This  divides  the  tree  into  connected  components.  Fig¬ 
ure  8  shows  how  the  program  aggregates  points  of  a  quasiperiodic  orbit  and 
recognizes  it  as  an  island  chain. 

To  compute  the  area  and  centroid  of  the  region  bounded  by  a  curve, 
KAM  generates  an  ordered  sequence  of  points  from  the  spanning  tree,  and 
spline-interpolates  the  sequence  to  obtain  a  smooth  curve.  Straightforward 
algorithms  are  then  applied  to  compute  the  area  and  centroid.  Shape  recog¬ 
nition  is  accomplished  using  scale- space  methods  pioneered  by  Witkin  [28]. 

2.2  AI  techniques  can  implement  deep  mathematical 
knowledge 

Viewed  as  abstract  examples  of  AI  technology,  our  demonstration  programs 
are  hardly  novel.  The  uniqueness  of  these  programs,  and  the  source  of  their 
power,  is  that  they  use  classic  AI  methods  to  exploit  specific  domain  knowl¬ 
edge  based  on  rigorous  mathematical  results. 

Plr  for  instance,  combines  geometric  reasoning,  symbolic  algebra,  and 
inequality  reasoning  to  test  whether  trajectories  of  a  piecewise  linear  system 
cross  between  adjacent  regions  in  phase  space.  For  a  trajectory  to  cross  from 
region  R  to  S  via  boundary  u,  its  tangent  t  at  the  intersection  point  with 
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Figure  7 :  Starting  with  the  successive  iterates  of  a  point,  KAM  classifies  orbits  using  al¬ 
gorithms  from  machine  vision.  As  shown  above,  a  quasiperiodic  orbit  can  be  distinguished 
from  a  chaotic  orbit  by  examining  the  branching  factor  in  the  Euclidean  minimal  spanning 
tree. 


Figure  8;  KAM  uses  the  minimal  spanning  tree  to  cluster  orbits  into  components.  The 
components  of  an  island  chain  can  be  isolated  by  detecting  long  edges  in  the  spanning  tree 
and  deleting  these  from  the  graph. 


u  must  form  an  acute  angle  with  the  normal  n,  eis  shown  in  Figure  9.  This 
geometric  condition  is  equivalent  to  the  algebraic  condition  that  the  inner 
product  t  n  be  positive.  Hence,  a  transition  exists  from  Rto  S  unless  t-n  <  0 
everywhere  on  u.  Plr  resolves  the  inequality  1  -  n  <  0  on  u  with  the  BOUNDED 
inequality  reasoner  (20j. 

Plr  combines  symbolic  reasoning  with  deep  knowledge  about  dynamical 
systems  to  interpret  the  transition  graphs  that  it  constructs.  For  example, 
the  transition  graph  for  the  van  der  Pol  equation  shows  that  all  trajectories 
spiral  eu-ound  the  origin,  but  tells  nothing  of  whether  they  move  inward,  move 
outward,  or  wobble  around.  Plr  invokes  a  difficult  theorem  to  prove  that 
all  trajectories  converge  to  a  unique  limit  cycle.  It  tests  the  preconditions  of 
the  theorem  by  proving  inequalities  and  manipulating  symbolic  expressions. 

The  KAM  program  limits  the  number  of  phase-space  trajectories  it  must 
explore  by  drawing  upon  constraint  analysis,  as  pioneered  by  Waltz  [27].  As 
in  any  constraint  analysis,  KAM  relies  upon  “grammar”  that  expresses  the 
consistent  ways  in  which  primitive  elements  can  be  combined.  In  KAM’s  case, 
the  primitive  elements  incorporate  sophisticated  mathematical  invariants. 
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Figure  9:  Trajectory  crossing  from  R  to  S  via  u:  the  tangent  t  at  the  crossing 
point  must  form  an  acute  angle  with  n,  the  normal  to  u  that  points  into  S. 

and  the  grammatical  rules  embody  deep  theorems  about  the  behavior  of 
dynamical  systems. 

One  such  invariant,  for  example,  is  the  rotation  number  of  an  orbit,  a 
quantity  that  measures  the  asymptotic  average  of  the  angular  distances  be¬ 
tween  any  two  successive  iterates  in  units  of  2ir-radians.  A  rule  in  RAM’s 
grammar  embodies  a  theorem  that  the  rotation  numbers  of  nearby  orbits 
must  change  continuously  [16].  As  an  example  of  how  this  is  used,  suppose 
that  RAM  has  located  two  nearby  almost-periodic  orbits  having  rotation 
numbers  pi  and  p2  respectively.  Suppose  pi  is  slightly  smaller  than  1/5,  and 
P2  slightly  larger.  With  only  these  two  orbits,  RAM’s  evolving  phase-space 
picture  cannot  be  complete.  By  continuity,  RAM  expects  to  find  a  third, 
nearby  orbit  with  rotation  number  exactly  equal  to  1  /5,  that  is,  a  periodic 
orbit  of  period  5,  which  RAM  proceeds  to  search  for  and  classify. 

In  a  similar  manner,  Abelson’s  Bifurcation  Interpreter  draws  upon  knowl¬ 
edge  of  the  geometry  of  typical  changes  in  the  steady-state  orbits  of  one- 
parameter  families  of  dynamical  systems.  Periodic  orbits  of  a  periodically- 
driven  oscillator  can  be  identified  as  fixed-points  of  the  period  map,  which 
maps  a  state  to  the  end-point  of  the  trajectory  starting  from  that  state  and 
evolving  for  one  period  of  the  drive.  Stability  of  an  orbit  is  determined  by 
the  stability  of  the  corresponding  fixed  point.  If  the  interpreter  notices  that 
a  stable  orbit  suddenly  becomes  unstable  as  the  family  parameter  increases, 
it  attempts  to  explain  this  change  as  the  result  of  a  bifurcation.  The  type 
of  bifurcation  can  be  conjectured  by  examining  the  eigenvalues  of  the  period 
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map  at  the  fixed  point,  and  the  conjecture  can  be  verified  by  a  seajch  that 
is  tailored  to  the  local  geometry  for  that  bifurcation  type. 

For  example,  for  a  stable  orbit,  the  complex  eigenvalues  of  the  correspond¬ 
ing  fixed  point  of  the  period  map  must  lie  within  the  unit  circle.  If  stability 
is  lost  with  an  eigenvalue  apparently  crossing  the  unit  circle  at  —1,  a  classi¬ 
fication  theorem  for  bifurcations  [11]  tells  the  interpreter  to  expect  that  this 
is  a  supercritical-flip  bifurcation,  which  corresponds  to  a  period  doubling  of 
the  orbit.  Near  the  bifurcation  one  should  expect  to  see  a  stable  orbit  just 
before  the  critical  parameter  value  and,  just  after  the  critical  value,  an  un¬ 
stable  orbit  together  with  a  stable  orbit  of  double  the  period.  For  this  type  of 
bifurcation,  the  interpreter  attempts  to  locate  the  new  expected  orbits  using 
a  search  technique  that  detects  fixed  points  by  computing  the  Poincare  index 
of  the  period  map  [12].  If  these  orbits  are  located,  the  bifurcation  is  probably 
a  flip.  If  a  different  local  geometry  is  found,  the  apparent  bifurcation  may  be 
the  result  of  numerical  error,  or  the  interaction  of  two  nearby  bifurcations, 
or  it  may  be  a  bifurcation  of  non-standard  type.  In  any  case,  the  result  of 
the  search  is  passed  to  a  critic  that  attempts  to  reconcile  the  local  results  for 
all  bifurcations  detected  and  produce  a  consistent  description. 

Going  beyond  general  knowledge  of  dynamical  systems,  the  Kineticist’s 
Workbench  program  employs  a  number  of  techniques  specific  to  the  domeiin 
of  chemical  kinetics.  For  example,  the  program  examines  the  qualitative  his¬ 
tory  of  a  reaction  simulation,  attempting  to  find  periods  of  time  during  which 
concentrations  of  some  species  may  be  treated  as  constant;  this  is  an  automa¬ 
tion  of  the  kind  of  steady-state  analysis  that  is  a  staple  of  kinetic  investiga¬ 
tion  [15].  Another  portion  of  the  program — that  portion  devoted  to  spotting 
fast  equilibria  within  a  mechanism — makes  extensive  use  of  the  reaction  net¬ 
work  formalism  developed  by  Feinberg,  Horn,  Jackson,  and  coworkers  at  the 
University  of  Rochester  [6].  An  especially  fruitful  result  of  their  work  is  the 
zero- deficiency  theorem,  which  provides  a  simple  algorithmic  test  for  deter¬ 
mining  whether  a  reaction  mechanism  gives  rise  to  stable  equilibria.  Finally, 
the  portion  of  the  Workbench  program  devoted  to  decomposing  mechanisms 
according  to  mutual  influence  between  species  may  be  used  to  identify  tran¬ 
sient  species,  and  to  drop  those  species  from  a  particular  simulation  as  soon 
as  their  concentrations  are  deemed  too  low  to  affect  the  remainder  of  the 
simulation. 


24 


(delinc-natvork  drivaa-van-dcr-pol 
((a  paraaatar  T/i‘3) 

(b  paraaatar  raiistanca) 

(d  driaa  voltaga)} 

(nl  &2  n3) 

(parts 

(nl-raa  non-linaar-raaiator  (a-*-  a3)  (a-  gad) 

(vie  (laabda  (v  i) 

(»  V  (-  (*  a  i  i  i)  (a  b  i)))))) 

(1  iadactor  (aa  al)  (a-  a2)) 

(e  ci^eltor  (aa  a2)  (a-  a3)) 

(a  voltaga-soarea  (aa  al)  (a-  gad)  (straagth  d)))) 


Figure  10:  The  wiring  diagram  of  a  simple  nonlinear  circuit  is  described  by  means  of 
a  special-purpose  language  of  electrical  networks.  This  description  can  be  automatically 
compiled  into  numerical  procedures  that  evolve  the  state  of  the  system. 

2.3  Numerical  experiments  can  be  prepared  automat¬ 
ically 

Translating  high-level  specifications  into  high-quality  numerical  routines  can 
be  a  tedious  and  error-prone  process  whose  difficulty  limits  the  utility  of  even 
the  most  powerful  numerical  computers.  A  program  by  H.  Abelson  and  G.J. 
Sussman  draws  upon  a  spectrum  of  computational  tools — numerical  meth¬ 
ods,  symbolic  dgebra,  and  semantic  constraints  (such  as  dimensions) — to 
automate  the  preparation  and  execution  of  numerical  simulations  [2].  These 
tools  are  designed  so  that  combined  methods,  tailored  to  particular  problems, 
can  be  constructed  on  the  fly.  One  can  use  symbolic  algebra  to  automatically 
generate  numerical  procedures,  and  one  can  use  domain-specific  constraints 
to  guide  algebraic  derivations  and  to  avoid  complexity. 

Figure  10  shows  the  wiring  diagram  for  a  simple  nonlinear  circuit,  a  driven 
van  der  Pol  oscillator  consisting  of  voltage  source,  a  capacitor,  an  inductor, 
and  a  nonlinear  resistor  with  the  cube-law  characteristic  v  =  ai^  —  bi.  The 
figure  also  shows  a  description  of  this  wiring  diagram  in  a  language  formu¬ 
lated  for  describing  electrical  networks.  This  description  specifies  circuit’s 
parameters,  its  primitive  parts,  and  how  the  parts  are  interconnected. 

Given  this  description,  the  program  combines  models  of  the  primitive 
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elements  to  form  equations  that  are  then  algebraically  solved  to  produce  state 
equations  for  the  van  der  Pol  oscillator.  The  state  equations  are  compiled  into 
a  procedure  (the  system  derivative)  that  will  evolve  the  system  numerically 
when  combined  with  an  appropriate  numerical  integrator. 

These  operations  can  involve  a  nontrivial  amount  of  adgebraiic  manipu¬ 
lation.  Even  for  systems  that  are  specified  in  closed  form,  most  nonlinear 
systems  cannot  be  algebraically  solved  to  produce  explicit  state  equations.  In 
the  general  case,  the  program  recognizes  variables  that  cannot  be  eliminated 
from  the  state  equations  and  compiles  an  iterative  scheme  for  approximating 
these  variables.  This  requires  symbolic  differentiation  to  produce  a  Jacobian 
that  is  incorporated  into  a  Newton-Raphson  search  and  to  augment  the  sys¬ 
tem  state  so  that  it  will  evolve  good  starting  points  for  Newton’s  method  at 
each  step  of  the  integration. 

For  applications  such  as  the  Bifurcation  Interpreter,  one  must  also  compile 
numerical  routines  that  find  period  orbits  and  track  them  as  the  system 
parameters  vary.  Finding  and  tracking  periodic  orbits  rests  upon  determining 
the  sensitivity  of  trajectories  to  variations  in  initial  state  and  parameters. 
This  can  be  done  by  evolving  the  variational  system,  obtained  by  symbolic 
manipulation  of  the  state  equations,  and  by  evolving  various  derived  systems 
obtained  by  differentiating  the  state  equations  with  respect  to  parameters. 
Figure  11  shows  a  numerical  routine,  automatically  generated  from  the  circuit 
description  in  figure  10,  that  can  be  combined  with  a  numerical  integrator 
to  evolve  the  states  and  the  variational  states  of  the  driven  van  der  Pol 
oscillator. 

The  result  of  the  physical  modeling,  algebraic  manipulation,  and  fmcy 
compilation  shown  above  in  figure  11  is  a  higher-order  procedure — a  system- 
derivative  generator  for  the  dynamical  system  under  study.  The  generator 
takes  as  arguments  numerical  values  for  the  system  parameters  and  produces 
a  system-derivative  procedure,  which  takes  a  system  state  vector  as  argument 
and  produces  a  differential  state  (a  vector  that  when  multiplied  by  an  incre¬ 
ment  of  time  is  an  increment  of  state).  This  system- derivative  procedure 
is  passed  to  an  integration  driver  that  returns  a  procedure  which,  given  an 
initial  state,  evolves  the  system  numerically. 

Since  all  system  derivative  procedures  are  constructed  to  respect  the  same 
conventional  interfaces,  we  may  choose  from  a  variety  of  integration  meth¬ 
ods.  Moreover,  the  integration  methods  themselves  can  be  automatically 
constructed  from  a  library  of  procedures  that  can  be  used  as  interchangeable 
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(laabda  (c.c  l.I  d  b  a) 

(laabda  (*varstate*) 

(let  ((t  (vactor-rel  ^varatate*  0)) 

(v.c  (vector-ref  *veu:stata*  1)) 

(i.l  (vector-ref  ♦varstate*  2)) 

(v.c. del. v.c  (vector-ref  avaratate*  3)) 

(v.c. del. i.l  (vector-ref  evaratate*  4)) 

(i.l. del. v.c  (vector-ref  ovaratate*  5)) 

(i.l. dal. i.l  (vector-ref  evaratate*  6))) 

(let  ((g27  (*  a  i.l  i.l))) 

(let  ((g28  (•  -3  g27))) 

(vector  1 

(/  i.l  c.c) 

(/  (+  (*  -1  g27  i.l)  (*  -1  v.c)  (•  b  i.l)  (d  t)) 

1.1) 

(/  v.c. del. i.l  c.c) 

(/  (+  (♦  “1  v.c. del. v.c) 

(«  b  v.c. del. i.l) 

(*  g28  v.c. del. i.l)) 

1.1) 

(/  i.l. del. i.l  c.c) 

(/  (+  (♦  “1  i.l. del. v.c) 

(e  b  i.l. del. i.l) 

(e  g28  i.l. del. i.l)) 

1.1))))))) 


Figure  11:  This  numerical  procedure  is  the  augmented  system  derivative  generator  for 
evolving  variational  states  for  the  driven  van  der  Pol  oscillator.  The  procedure  was  auto¬ 
matically  generated  from  the  circuit  description  shown  above. 
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components  in  the  construction  of  traditional  applications.  Going  further, 
we  believe  that  it  is  not  difficult  to  automatically  implement  these  numerical 
procedures  as  special-purpose  hardware,  like  the  Orrery. 

To  support  the  automatic  construction  of  numerical  procedures,  we  are 
developing  a  kernel  numerical  library  that  is  organized  around  the  liberal 
use  of  high-order  procedural  abstractions.  For  example,  figure  12  illus¬ 
trates  this  mix-and-match  construction  of  numerical  routines,  expressing 
Romberg’s  method  of  quadrature  as  a  combination  of  trapezoidal  integration 
and  Richardson  extrapolation,  following  the  exposition  given  by  M.  Haifant 
and  G.J.  Sussman  in  [8].  Such  a  formulation  is  valuable  in  that  it  separates 
the  ideas  into  several  independent  pieces.  Clever  ideas  need  be  coded  and 
debugged  only  once,  in  a  context  independent  of  the  particular  application, 
thus  enhancing  the  reliability  of  software  built  in  this  way.  For  instance, 
G.  Roylance  [19]  shows  how  to  construct  high-performance  implementations 
of  special  functions,  abstracting  recurrent  themes  such  as  Chebyshev  econ¬ 
omization.  His  automatically-constructed  procedure  for  computing  Bessel 
functions  is  40  times  more  accurate,  for  the  same  number  of  terms,  than  the 
approximation  specified  in  the  National  Bureau  of  Standards  tables  [3].  More 
significantly,  Roylance’s  formulation  clearly  exposes  the  underlying  approx¬ 
imation  methods  so  that  parameters,  such  as  the  required  precision  of  the 
routines,  can  be  changed  at  will. 

Besides  providing  a  convenient  target  for  automatic  construction  of  nu¬ 
merical  procedures,  powerful  abstraction  mechanisms  help  us  to  express  some 
of  the  vocabulary  and  methods  of  numerical  analysis  in  a  form  that  is  close 
to  the  mathematical  theory,  and  is  thus  easy  to  understand  and  check.  A 
program  is  a  communication,  not  just  between  programmers  and  computers, 
but  also  between  programmers  and  human  readers  of  the  program;  quite  of¬ 
ten,  between  the  programmer  and  him/herself.  One  power  of  programming  is 
that  it  allows  one  to  make  the  knowledge  of  methods  explicit,  so  that  meth¬ 
ods  can  be  studied  as  theoretical  entities.  Traditional  numerical  programs 
are  hand-crafted  for  each  application.  The  traditional  style  does  not  admit 
such  explicit  decomposition  and  naming  of  methods,  thus  forfeiting  much  of 
the  power  and  joy  of  programming. 
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Figure  12:  Romberg’s  method  of  quadrature  can  be  built  by  combining  a  primitive  trape¬ 
zoidal  integrator  with  an  accelerator  that  speeds  convergence  of  sequences  by  Richardson 
extrapolation.  The  result  is  an  infinite  sequence  (stream)  of  increasingly  accurate  approx¬ 
imations  to  the  definite  integral.  The  same  Richardson  accelerator  can  be  combined  with 
other  sequence  generators  to  build  other  classical  numerical  routines. 
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3  Intelligent  tools  are  feasible 

The  work  described  in  the  preceding  sections  demonstrates  much  of  the  tech¬ 
nology  required  to  produce  programs  that  can  serve  scientists  and  engineers 
eis  intelligent  problem-solving  partners,  programs  such  as  the  engineer’s  as¬ 
sistant  that  we  envisioned  at  the  beginning  of  this  paper. 

We  have  shown  how  to  use  symbolic  algebra  to  compile  high-level  descrip¬ 
tions  such  as  circuit  diagrams  directly  into  numerical  modeling  and  simula¬ 
tion  programs  whose  elements  can  be  automatically  generated  from  a  librau'y 
of  mix-and-match  numerical  subroutines  expressed  at  appropriate  levels  of 
abstraction.  Our  experience  with  the  Digital  Orrery  proves  that  such  nu¬ 
merical  programs  can  be  run  at  supercomputer  speeds,  without  the  cost  of  a 
general-purpose  supercomputer.  The  Bifurcation  Interpreter,  KAM  and  Plr 
demonstrate  that  intelligent  programs  incorporating  knowledge  of  dynami¬ 
cal  systerhs  can  automatically  control  and  monitor  numerical  experiments 
<ind  interpret  the  results  in  qualitative  terms.  The  Kineticist’s  Workbench 
illustrates  how  these  capabilities  can  be  combined  with  knowledge  about  a 
particular  domain  to  produce  a  sophisticated  tool  for  modeling  and  analysis. 

3.1  Higher-dimensional  systems  are  hard 

Most  systems  of  interest  have  more  than  two  degrees  of  freedom,  yet  the 
KAM  and  Plr  programs  and  the  Bifurcation  Interpreter  depend  upon  spe¬ 
cial  properties  of  low-dimensional  systems.  The  granunar  of  possible  phase 
portraits  and  the  catalog  of  generic  bifurcations  embodied  in  these  progreims 
cannot  easily  be  extended  to  higher  dimensional  systems.  On  the  other  hand, 
there  are  qualitative  features  of  such  systems  that  can  be  usefully  extracted 
and  used  to  guide  numerical  experiments. 

Exploring  the  qualitative  behavior  of  high-dimensional  systems  requires  a 
combination  of  analytic  and  numeric  methods.  Analytic  methods  can  provide 
clear  definitive  information,  but  are  often  hard  to  apply  or  unavailable.  There 
are  well-established  method  for  deriving  the  local  behavior  of  trajectories  in 
the  neighborhood  of  fixed  points,  but  few  tools  exist  for  determining  global 
behavior.  On  the  other  hand,  a  program  could  detect  a  saddle  analytically, 
cadculate  its  stable  and  unstable  manifolds  numerically,  determine  whether 
the  manifolds  intersect  each  other,  and  draw  conclusions  about  global  behav¬ 
ior. 
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Moreover,  even  in  high-dimensional  spaces,  it  is  still  possible  to  use  clus¬ 
tering  techniques  to  examine  the  set  of  iterates  of  a  map  or  the  flow  of 
a  differential  equation  and  determine  if  a  trajectory  is  confined  to  a  lower 
dimensional  submanifold  of  the  formal  state  space.  Eax:h  reduction  in  dimen¬ 
sion  is  evidence  of  an  integral  of  motion,  such  as  conservation  of  energy.  One 
can  also  (as  people  do)  apply  visual  recognition  techniques  to  low- dimensional 
sections  and  projections  of  the  full  space.  Despite  the  fact  that  orbit  types 
and  bifurcations  in  high- dimensional  spaces  have  not  been  completely  classi¬ 
fied,  nonetheless  it  is  still  possible  to  recognize  qualitatively  different  regions 
of  behavior,  and  to  map  out  these  regions  in  state  space  and  parameter  space. 

3.2  Computers,  like  people,  need  imagistic  reasoning 

In  observing  professional  physicists  and  engineers,  we  are  often  struck  by 
how  an  expert’s  “intuitive  grasp”  of  a  field  is  hard  to  articulate  verbally. 
This  is  perhaps  indicative  of  the  use  of  non-verbal  reasoning  processes  as 
part  of  the  process  of  solving  otherwise  verbally  presented  problems.  We 
observe  scientists,  mathematicians,  and  engineers  continually  using  graphical 
representations  to  organize  their  thoughts  about  a  problem.  The  programs  we 
are  developing  use  numerical  methods  as  a  means  of  shifting  back  and  forth 
between  symbolic  and  geometric  methods  of  reasoning.  The  programs  not 
only  draw  graphs  and  state-space  diagrams,  but  they  look  at  these  diagrams 
and  hold  them  in  their  “mind’s  eye”  so  that  powerful  visual  mechanisms  can 
be  brought  to  bear  on  what  otherwise  would  be  purely  symbolic  problems. 

The  idea  that  problem  solvers  employing  visual,  analogue,  or  diagram¬ 
matic  representations  can  be  more  effective  than  those  relying  on  linguistic 
representations  alone  is  not  new.  Even  before  1960,  Gelemter’s  Geometry- 
Theorem  Proving  Machine  [10]  used  diagrams  to  filter  goals  generated  by 
backward  chaining.  Nevins’s  [18]  forward-chaining  theorem  prover  focussed 
its  forward  deduction  of  facts  on  those  lines  explicitly  drawn  in  a  diagram. 
St2dlman  and  Sussman’s  EL  [23]  program  performed  antecedent  deductions 
in  circuit  analysis  by  exploiting  the  finite  connectivity  of  devices. 

What  is  provocative,  however,  is  the  suggestion  that  our  thought  pro¬ 
cesses  are  importantly  imagistic,  and  that  visual  thinking  may  play  a  crucial 
role  in  problem  solving.  In  scientific  computation  there  has  been  tremendous 
emphasis  on  visualization,  but  this  has  mostly  meant  the  development  of 
computer-graphics  technology  to  aid  human  visualization  [17].  We  believe 
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that  imagistic  reasoning  is  a  very  general  class  of  problem-solving  strategies, 
each  with  its  own  appropriate  representations  and  technical  support.  The 
programs  discussed  above  suggest  that  it  may  be  at  least  as  important  for 
scientific  computation  to  develop  visualization  aids  for  programs  as  well  as 
for  people. 
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